/*
 * Parallel Force-Directed Edge Bundling using CUDA.
 *
 * Author: jasmine2800@hotmail.com
 *
 */


#include <stdio.h>
#include <string>
#include <math.h>
#include "time.h"

#define BLOCK_WIDTH 16
#define BLOCK_SIZE_1 32
#define numEdges 465
#define numEdgesEx 480
#define numSubdivs 11
#define numPoints 31

#define EPS  1e-7

#define Macro_updateViewAfterEachStep 1
#define Macro_numCycles 10
#define Macro_P 1
#define Macro_S 1.5
#define Macro_I 150
#define Macro_K 0.1
#define Macro_repulsionAmount 1.0
#define Macro_stepDampingFactor 0.5
#define Macro_edgeCompatibilityThreshold 0.6
#define Macro_directionAffectsCompatibility 1
#define Macro_binaryCompatibility 0
#define Macro_useInverseQuadraticModel 0
#define Macro_useRepulsionForOppositeEdges 1
#define Macro_useSimpleCompatibilityMeasure 0
#define Macro_edgeValueAffectsAttraction /*1*/0
#define Macro_subdivisionPointsCycleIncreaseRate 1.3



//********************************************************************
typedef struct EP
{
    float x, y;
} Point;

float *edgeLengths_h;
float *edgeLengths_d;

float *edgeValues_h;
float *edgeValues_d;

float edgeValueMax, edgeValueMin;

Point endPoints[numPoints];

float *data_x_h;
float *data_y_h;
float *data_x_d;
float *data_y_d;

float *compatibility_h;
float *compatibility_d;

float *Fsi_x_d,*Fsi_y_d,//to store position deviation caused by spring force.
      *Fei_x_d,*Fei_y_d,//to store position deviation caused by electrostatic force.
      *Fei_x_d_0,*Fei_y_d_0;//the intermediate result of Fei_d. //Fei_d=sum(Fei_d_0).

int P;    // number of subdivision points (will increase with every cycle).
float Pfloat; // used to keep the float value to keep the stable increase rate.
float S;   // step size.
int I;    // number of iteration steps performed during a cycle.



//********************************************************************
#define computeDistance(x1,y1,x2,y2) sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))



//********************************************************************
__device__ float signum(float x)
{
    if (x == 0.0) return 0;
    else if (x < 0.0) return -1.0;
    else return 1.0;
}



//********************************************************************
//compute spring force.
__global__ void computeSpringForce(float S,int P,float *edgeLengths_d,float *data_x_d,float *data_y_d,float *Fsi_x_d,float *Fsi_y_d)
{
    __shared__ float data_x_s[numSubdivs+2][BLOCK_WIDTH];
    __shared__ float data_y_s[numSubdivs+2][BLOCK_WIDTH];

    __shared__ float Fsi_x_s[numSubdivs+1][BLOCK_WIDTH];
    __shared__ float Fsi_y_s[numSubdivs+1][BLOCK_WIDTH];

    __shared__ float k_p_s[BLOCK_WIDTH];//to store the 16 edge's spring rates, used by current block.

    int yIndex=16*blockIdx.y+threadIdx.x;

    if(threadIdx.y==0)   //the 1st half warp do this.
    {
        k_p_s[threadIdx.x] = Macro_K / (edgeLengths_d[16*blockIdx.y+threadIdx.x] * (P+1));
        //read end points.
        data_x_s[P+1][threadIdx.x]=data_x_d[blockIdx.y*16+threadIdx.x];
        data_y_s[P+1][threadIdx.x]=data_y_d[blockIdx.y*16+threadIdx.x];
    }

    if(/*yIndex<numEdges&&*/threadIdx.y<P+1)//each line of threads read one column of points.
    {
        //read points.
        data_x_s[threadIdx.y][threadIdx.x]=data_x_d[(threadIdx.y+1)*numEdgesEx+blockIdx.y*16+threadIdx.x];
        data_y_s[threadIdx.y][threadIdx.x]=data_y_d[(threadIdx.y+1)*numEdgesEx+blockIdx.y*16+threadIdx.x];
    }
    __syncthreads();

    if(/*yIndex<numEdges&&*/threadIdx.y<P+1)//each line of threads compute one column of spring force.
    {
        Fsi_x_s[threadIdx.y][threadIdx.x]=data_x_s[threadIdx.y+1][threadIdx.x]-data_x_s[threadIdx.y][threadIdx.x];
        Fsi_y_s[threadIdx.y][threadIdx.x]=data_y_s[threadIdx.y+1][threadIdx.x]-data_y_s[threadIdx.y][threadIdx.x];

        if (abs(k_p_s[threadIdx.x]) < 1.0)
        {
            Fsi_x_s[threadIdx.y][threadIdx.x] *= k_p_s[threadIdx.x];
            Fsi_y_s[threadIdx.y][threadIdx.x] *= k_p_s[threadIdx.x];
        }
    }
    __syncthreads();

    //if(yIndex<numEdges)
    {
        //sum up to yield the spring force taken by each point. 
        if (threadIdx.y>0&&threadIdx.y<P+1)
        {
            Fsi_x_d[(threadIdx.y-1)*numEdgesEx+blockIdx.y*16+threadIdx.x]=data_x_s[threadIdx.y][threadIdx.x]+
                    (Fsi_x_s[threadIdx.y][threadIdx.x]-
                     Fsi_x_s[threadIdx.y-1][threadIdx.x])*S;

            Fsi_y_d[(threadIdx.y-1)*numEdgesEx+blockIdx.y*16+threadIdx.x]=data_y_s[threadIdx.y][threadIdx.x]+
                    (Fsi_y_s[threadIdx.y][threadIdx.x]-
                     Fsi_y_s[threadIdx.y-1][threadIdx.x])*S;
        }
    }
}



//********************************************************************
//sum up the columns of the matrix of electrostatic force, yield the electrostatic force taken by each point.
__global__ void sumElectroForce(
    int P,
    float *Fei_x_d_0,
    float *Fei_y_d_0,
    float *Fei_x_d,
    float *Fei_y_d)
{
    __shared__ float Fei_x_s[16][17];//electrostatic force.
    __shared__ float Fei_y_s[16][17];

    __shared__ float sum_x_s[16];
    __shared__ float sum_y_s[16];

    int xIndex=blockDim.x*blockIdx.x+threadIdx.x;
    int yIndex=blockDim.y*blockIdx.y+threadIdx.y;
    int counts;
    int s;
    int i;

    for (i=0; i<P; i++)//iterate to handle the columns one by one.
    {
        if (yIndex<numEdges)
        {
            if (threadIdx.x==0)
            {
                sum_x_s[threadIdx.y]=0;
                sum_y_s[threadIdx.y]=0;
            }
            for (counts=480/16; counts>0; counts-=16,xIndex+=16)
            {
                //load data into share mem.
                //if(counts<16)
                {
                    Fei_x_s[threadIdx.y][threadIdx.x]=0;
                    Fei_y_s[threadIdx.y][threadIdx.x]=0;
                }
                if(threadIdx.x<counts)
                {
                    Fei_x_s[threadIdx.y][threadIdx.x]=Fei_x_d_0[i*numEdgesEx*numEdgesEx+yIndex*numEdgesEx+xIndex];
                    Fei_y_s[threadIdx.y][threadIdx.x]=Fei_y_d_0[i*numEdgesEx*numEdgesEx+yIndex*numEdgesEx+xIndex];
                }
                //sum up.
                for (s=8; s>0; s>>=1)
                {
                    if (threadIdx.x<s)
                    {
                        Fei_x_s[threadIdx.y][threadIdx.x]+=Fei_x_s[threadIdx.y][threadIdx.x+s];
                        Fei_y_s[threadIdx.y][threadIdx.x]+=Fei_y_s[threadIdx.y][threadIdx.x+s];
                    }
                }
                if (threadIdx.x==0)
                {
                    sum_x_s[threadIdx.y]+=Fei_x_s[threadIdx.y][0];
                    sum_y_s[threadIdx.y]+=Fei_y_s[threadIdx.y][0];
                }
            }
        }
        __syncthreads();//barrier.

        if (yIndex<numEdges)
        {
            //write back to global mem.
            if (threadIdx.y==0)
            {
                Fei_x_d[i*numEdgesEx+blockDim.y*blockIdx.y+threadIdx.x]=sum_x_s[threadIdx.x];
                Fei_y_d[i*numEdgesEx+blockDim.y*blockIdx.y+threadIdx.x]=sum_y_s[threadIdx.x];
            }
        }
        __syncthreads();//barrier.
    }
}



//********************************************************************
//compute intermediate result of electrostatic force.
__global__ void computeElectroForce(
    int P,
    float S,
    float *compatibility_d,
    float *data_x_d,float *data_y_d,
    float *Fei_x_d_0,float *Fei_y_d_0,
    float* edgeValues_d,
    float edgeValueMax,float edgeValueMin)
{
    __shared__ float cd_s[16][16];//compatibility.
    __shared__ float edgeValues_s[2][16];
    __shared__ float data_x_s[2][16];//matrix of points.
    __shared__ float data_y_s[2][16];
    __shared__ float Fei_x_s[16][17];//to store electrostatic force (16*16), the extra column is to avoid bank conflict.
    __shared__ float Fei_y_s[16][17];


    int xIndex=blockDim.x*blockIdx.x+threadIdx.x;
    int yIndex=blockDim.y*blockIdx.y+threadIdx.y;
    int i,j;

    float sum_x;
    float sum_y;
    float v_x;
    float v_y;
    float m;
    float d;

    if (blockIdx.x<=blockIdx.y)   //only the blocks under the diagonal do this.
    {
        //load compatibility.
        cd_s[threadIdx.y][threadIdx.x]=compatibility_d[yIndex*numEdgesEx+xIndex];

        if (threadIdx.y==0) //the 1st warp does this.
        {
            edgeValues_s[0][threadIdx.x]=edgeValues_d[xIndex];
        }
        else if (threadIdx.y==2) //the 2nd warp does this.
        {
            edgeValues_s[1][threadIdx.x]=edgeValues_d[blockDim.y*blockIdx.y+threadIdx.x];
        }
        __syncthreads();//barrier.

        for (i=0; i<P; i++)   //iterate to handle the columns one by one.
        {
            //load data into share mem.
            if (threadIdx.y==0)   //the 1st warp does this.
            {
                data_x_s[0][threadIdx.x]=data_x_d[i*numEdgesEx+xIndex];
                data_y_s[0][threadIdx.x]=data_y_d[i*numEdgesEx+xIndex];
            }
            else if (threadIdx.y==2)     //the 2nd warp does this.
            {
                data_x_s[1][threadIdx.x]=data_x_d[i*numEdgesEx+blockDim.y*blockIdx.y+threadIdx.x];
                data_y_s[1][threadIdx.x]=data_y_d[i*numEdgesEx+blockDim.y*blockIdx.y+threadIdx.x];
            }

            Fei_x_s[threadIdx.y][threadIdx.x]=0;
            Fei_y_s[threadIdx.y][threadIdx.x]=0;

            __syncthreads();//barrier.

            //begin to compute electrostatic force.
            if (xIndex<yIndex)
            {
                if(abs(cd_s[threadIdx.y][threadIdx.x])>=Macro_edgeCompatibilityThreshold)
                {
                    v_x=data_x_s[0][threadIdx.x]-data_x_s[1][threadIdx.y];
                    v_y=data_y_s[0][threadIdx.x]-data_y_s[1][threadIdx.y];
                    d=sqrt(v_x*v_x+v_y*v_y);
                    if (Macro_useInverseQuadraticModel)
                    {
                        m=(cd_s[threadIdx.y][threadIdx.x]/d)/(d*d);
                    }
                    else
                    {
                        m=(cd_s[threadIdx.y][threadIdx.x]/d)/d;
                    }
                    if(cd_s[threadIdx.y][threadIdx.x]<0)
                    {
                        m*=Macro_repulsionAmount;
                    }
                    if(Macro_edgeValueAffectsAttraction)
                    {
                        m*=(1.0 + max(-1.0, (edgeValues_s[0][threadIdx.x] - edgeValues_s[1][threadIdx.y])/(edgeValueMax + edgeValueMin)));
                    }
                    if(abs(m*S)>1.0)
                    {
                        m=signum(m)/S;
                    }
                    v_x*=m;
                    v_y*=m;
                    Fei_x_s[threadIdx.y][threadIdx.x]=v_x;
                    Fei_y_s[threadIdx.y][threadIdx.x]=v_y;
                    if(blockIdx.x==blockIdx.y)//the block on diagonal.
                    {
                        Fei_x_s[threadIdx.x][threadIdx.y]=-v_x;
                        Fei_y_s[threadIdx.x][threadIdx.y]=-v_y;
                    }
                }
            }

            __syncthreads();//barrier.

            //sum up the columns.
            if (threadIdx.y==0)   //the 1st warp does this.
            {
                sum_x=0;
                sum_y=0;
                for (j=0; j<16; j++)
                {
                    sum_x+=Fei_x_s[threadIdx.x][j];
                    sum_y+=Fei_y_s[threadIdx.x][j];
                }
                //write back to global mem.
                Fei_x_d_0[i*numEdgesEx*numEdgesEx+(blockDim.y*blockIdx.y+threadIdx.x)*numEdgesEx+blockIdx.x]=sum_x*S;
                Fei_y_d_0[i*numEdgesEx*numEdgesEx+(blockDim.y*blockIdx.y+threadIdx.x)*numEdgesEx+blockIdx.x]=sum_y*S;
            }
			//sum up the rows.
            else if (blockIdx.x<blockIdx.y&&threadIdx.y==2) //the 2nd warp does this.
            {
                sum_x=0;
                sum_y=0;
                for (j=0; j<16; j++)
                {
                    sum_x+=Fei_x_s[j][threadIdx.x];
                    sum_y+=Fei_y_s[j][threadIdx.x];
                }
                //write back to global mem.
                Fei_x_d_0[i*numEdgesEx*numEdgesEx+(blockDim.x*blockIdx.x+threadIdx.x)*numEdgesEx+blockIdx.y]=-sum_x*S;
                Fei_y_d_0[i*numEdgesEx*numEdgesEx+(blockDim.x*blockIdx.x+threadIdx.x)*numEdgesEx+blockIdx.y]=-sum_y*S;
            }
        }
    }
}



//********************************************************************
//Fsi stores the sum up of original point and spring force.
//Fei stores the electrostatic force.
//so, just sum up Fei and Fei to get the updated point.
__global__ void updateSubdivs(int P,float *data_x_d,float *data_y_d,float *Fsi_x_d,float *Fsi_y_d,float *Fei_x_d,float *Fei_y_d)
{
    int yIndex=blockDim.y*blockIdx.y+threadIdx.y;

    if (yIndex<P*numEdgesEx)
    {
        data_x_d[yIndex]=Fsi_x_d[yIndex]+Fei_x_d[yIndex];
        data_y_d[yIndex]=Fsi_y_d[yIndex]+Fei_y_d[yIndex];
    }

}



//********************************************************************
//compute edge compatibility.
__global__ void calcEdgeCompatibility(float *data_x_d,float *data_y_d,float *edgeLengths_d,float *compatibility_d)
{

    //the coordinate of current thread.
    int xIndex=blockDim.x*blockIdx.x+threadIdx.x;
    int yIndex=blockDim.y*blockIdx.y+threadIdx.y;
    int tIndex=blockDim.y*blockIdx.y+threadIdx.x;

    __shared__ float edgeStarts_x_s[2][BLOCK_WIDTH];
    __shared__ float edgeStarts_y_s[2][BLOCK_WIDTH];

    __shared__ float edgeEnds_x_s[2][BLOCK_WIDTH];
    __shared__ float edgeEnds_y_s[2][BLOCK_WIDTH];

    __shared__ float edgeLengths_s[2][BLOCK_WIDTH];

    if(blockIdx.x<=blockIdx.y)
    {
        //load the start point into share mem.
        if(blockIdx.x==blockIdx.y)//only the blocks on the diagonal do this.
        {
            if (threadIdx.y==0)   //the 1st warp do this.
            {
                edgeStarts_x_s[0][threadIdx.x]=data_x_d[numEdgesEx+xIndex];
                edgeStarts_y_s[0][threadIdx.x]=data_y_d[numEdgesEx+xIndex];

                edgeEnds_x_s[0][threadIdx.x]=data_x_d[xIndex];
                edgeEnds_y_s[0][threadIdx.x]=data_y_d[xIndex];

                edgeLengths_s[0][threadIdx.x]=edgeLengths_d[xIndex];

                edgeStarts_x_s[1][threadIdx.x]=edgeStarts_x_s[0][threadIdx.x];
                edgeStarts_y_s[1][threadIdx.x]=edgeStarts_y_s[0][threadIdx.x];

                edgeEnds_x_s[1][threadIdx.x]=edgeEnds_x_s[0][threadIdx.x];
                edgeEnds_y_s[1][threadIdx.x]=edgeEnds_y_s[0][threadIdx.x];

                edgeLengths_s[1][threadIdx.x]=edgeLengths_s[0][threadIdx.x];
            }
        }
        else//the blocks under the diagonal do this.
        {
            if (threadIdx.y==0)   //the 1st warp do this.
            {

                edgeStarts_x_s[0][threadIdx.x]=data_x_d[numEdgesEx+xIndex];
                edgeStarts_y_s[0][threadIdx.x]=data_y_d[numEdgesEx+xIndex];

                edgeEnds_x_s[0][threadIdx.x]=data_x_d[xIndex];
                edgeEnds_y_s[0][threadIdx.x]=data_y_d[xIndex];

                edgeLengths_s[0][threadIdx.x]=edgeLengths_d[xIndex];

            }
            else if(threadIdx.y==2) //the 2nd warp do this.
            {
                edgeStarts_x_s[1][threadIdx.x]=data_x_d[numEdgesEx+tIndex];
                edgeStarts_y_s[1][threadIdx.x]=data_y_d[numEdgesEx+tIndex];

                edgeEnds_x_s[1][threadIdx.x]=data_x_d[tIndex];
                edgeEnds_y_s[1][threadIdx.x]=data_y_d[tIndex];

                edgeLengths_s[1][threadIdx.x]=edgeLengths_d[tIndex];
            }
        }

        __syncthreads();//barrier.

        if(yIndex<numEdges&&xIndex<yIndex)
        {
            float C;
            float Cv;

            float x1,x2,y1,y2,x,y,x4,y4;
            float tmp;
			
            x1=edgeEnds_x_s[0][threadIdx.x]-edgeStarts_x_s[0][threadIdx.x];
            y1=edgeEnds_y_s[0][threadIdx.x]-edgeStarts_y_s[0][threadIdx.x];
            x2=edgeEnds_x_s[1][threadIdx.y]-edgeStarts_x_s[1][threadIdx.y];
            y2=edgeEnds_y_s[1][threadIdx.y]-edgeStarts_y_s[1][threadIdx.y];
            tmp=edgeLengths_s[0][threadIdx.x];
            tmp*=edgeLengths_s[1][threadIdx.y];
			
            // angle compatibility.
            if(Macro_directionAffectsCompatibility)
            {
                C = ((x1*x2+y1*y2) / tmp + 1.0) / 2.0;
            }
            else
            {
                C = abs((x1*x2+y1*y2) / tmp);
            }
            if (abs(C) < EPS)
            {
                C = 0.0;    // this led to errors (when Ca == -1e-12).
            }
            if (abs(abs(C) - 1.0) < EPS)
            {
                C = 1.0;
            }
            if (Macro_useRepulsionForOppositeEdges)
            {
                tmp = (x1*x2+y1*y2) /(sqrt(x1 * x1 + y1 * y1) * sqrt(x2 * x2 + y2 * y2));
                if (tmp < 0)
                {
                    C = -C;
                }
            }

            // scale compatibility.
            tmp = (edgeLengths_s[0][threadIdx.x] + edgeLengths_s[1][threadIdx.y])/2;
            C *=( 2 / (
                      (tmp / min(edgeLengths_s[0][threadIdx.x], edgeLengths_s[1][threadIdx.y]))  +
                      (max(edgeLengths_s[0][threadIdx.x], edgeLengths_s[1][threadIdx.y]) / tmp)
                  ));

            // position compatibility.
            x1=(edgeEnds_x_s[0][threadIdx.x] + edgeStarts_x_s[0][threadIdx.x])/2;
            y1=(edgeEnds_y_s[0][threadIdx.x] + edgeStarts_y_s[0][threadIdx.x])/2;
            x2=(edgeEnds_x_s[1][threadIdx.y] + edgeStarts_x_s[1][threadIdx.y])/2;
            y2=(edgeEnds_y_s[1][threadIdx.y] + edgeStarts_y_s[1][threadIdx.y])/2;
            C *= (tmp / (tmp + sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))));

            // visibility compatibility.
            if (abs(C) >0.9)
            {
                x1=edgeStarts_x_s[0][threadIdx.x];
                y1=edgeStarts_y_s[0][threadIdx.x];
                x2=edgeEnds_x_s[0][threadIdx.x];
                y2=edgeEnds_y_s[0][threadIdx.x];
                x=edgeStarts_x_s[1][threadIdx.y];
                y=edgeStarts_y_s[1][threadIdx.y];
                tmp = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
                //i0
                x4=x1 + (x2-x1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);
                y4=y1 + (y2-y1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);

                x=edgeEnds_x_s[1][threadIdx.y];
                y=edgeEnds_y_s[1][threadIdx.y];
                tmp = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
                //i1
                x=x1 + (x2-x1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);
                y=y1 + (y2-y1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);

                x1=(x + x4)/2;
                y1=(y + y4)/2;

                x2=(edgeEnds_x_s[0][threadIdx.x] + edgeStarts_x_s[0][threadIdx.x])/2;
                y2=(edgeEnds_y_s[0][threadIdx.x] + edgeStarts_y_s[0][threadIdx.x])/2;

                Cv=1 - 2 * sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) / sqrt((x-x4)*(x-x4) + (y-y4)*(y-y4));
                if(Cv<0)
                {
                    Cv=0;
                }
				
                x1=edgeStarts_x_s[1][threadIdx.y];
                y1=edgeStarts_y_s[1][threadIdx.y];
                x2=edgeEnds_x_s[1][threadIdx.y];
                y2=edgeEnds_y_s[1][threadIdx.y];
                x=edgeStarts_x_s[0][threadIdx.x];
                y=edgeStarts_y_s[0][threadIdx.x];
                tmp = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
                //i0
                x4=x1 + (x2-x1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);
                y4=y1 + (y2-y1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);

                x=edgeEnds_x_s[0][threadIdx.x];
                y=edgeEnds_y_s[0][threadIdx.x];
                tmp = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
                //i1
                x=x1 + (x2-x1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);
                y=y1 + (y2-y1)*((y1-y)*(y1-y2) - (x1-x)*(x2-x1)) / (tmp * tmp);

                x1=(x + x4)/2;
                y1=(y + y4)/2;
                x2=(edgeEnds_x_s[1][threadIdx.y] + edgeStarts_x_s[1][threadIdx.y])/2;
                y2=(edgeEnds_y_s[1][threadIdx.y] + edgeStarts_y_s[1][threadIdx.y])/2;

                tmp=1 - 2 * sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) / sqrt((x-x4)*(x-x4) + (y-y4)*(y-y4));
                if(tmp<0)
                {
                    tmp=0;
                }
                if(tmp<Cv)
                {
                    Cv=tmp;
                }
            }
            else
            {
                Cv = 1.0;
            }

            C*= Cv;
            compatibility_d[numEdgesEx*BLOCK_WIDTH*blockIdx.y+BLOCK_WIDTH*blockIdx.x+numEdgesEx*threadIdx.y+threadIdx.x]=C;
        }
    }
}



//********************************************************************
//make preparation for calcEdgeCompatibility.
void calcEdgeCompatibilityMeasures()
{

    dim3 dimBlock(BLOCK_WIDTH,BLOCK_WIDTH);
    dim3 dimGrid(numEdgesEx/BLOCK_WIDTH,numEdgesEx/BLOCK_WIDTH);

    //allocate CompatibleEdge on the device.
    size_t size=numEdgesEx*numEdgesEx*sizeof(float);
    cudaMalloc((void**)&compatibility_d,size);

    compatibility_h=(float*)malloc(size);

    calcEdgeCompatibility<<<dimGrid,dimBlock>>>(data_x_d,data_y_d,edgeLengths_d,compatibility_d);
    cudaThreadSynchronize();

    cudaMemcpy(compatibility_h,compatibility_d,size,cudaMemcpyDeviceToHost);
}



//********************************************************************
//add subdivision points.
//if cycle==0.
__global__ void addMidPoints(float *data_x_d,float *data_y_d)
{
    int yIndex=blockDim.y*blockIdx.y+threadIdx.y;

    data_x_d[2*numEdgesEx+yIndex]=(data_x_d[yIndex]+data_x_d[numEdgesEx+yIndex])*0.5;
    data_y_d[2*numEdgesEx+yIndex]=(data_y_d[yIndex]+data_y_d[numEdgesEx+yIndex])*0.5;
}



//********************************************************************
//add subdivision points.
//if cycle>0.
__global__ void addSubdivisionPoints(int P,int prevP,float *data_x_d,float *data_y_d)
{
    int tIndex=16*blockIdx.y+threadIdx.x;

    __shared__ float data_x_s[numSubdivs+2][BLOCK_WIDTH];
    __shared__ float data_y_s[numSubdivs+2][BLOCK_WIDTH];

    // bigger array for subdivision points of the next cycle.
    __shared__ float newEdgePoints_x[numSubdivs][BLOCK_WIDTH];
    __shared__ float newEdgePoints_y[numSubdivs][BLOCK_WIDTH];

    __shared__ float segmentLen[numSubdivs+1][BLOCK_WIDTH];

    //load data into share mem.
    if(threadIdx.y==0)   //the 1st half warp do this.
    {
        //read end points.
        data_x_s[prevP+1][threadIdx.x]=data_x_d[tIndex];
        data_y_s[prevP+1][threadIdx.x]=data_y_d[tIndex];
    }

    if(tIndex<numEdges&&threadIdx.y<=prevP)//each line of threads read one column of points.
    {
        data_x_s[threadIdx.y][threadIdx.x]=data_x_d[(threadIdx.y+1)*numEdgesEx+tIndex];
        data_y_s[threadIdx.y][threadIdx.x]=data_y_d[(threadIdx.y+1)*numEdgesEx+tIndex];
    }
    __syncthreads();

    if(tIndex<numEdges&&threadIdx.y<=prevP)//each line of threads compute one column of lengths.
    {
        segmentLen[threadIdx.y][threadIdx.x]=computeDistance(
                data_x_s[threadIdx.y][threadIdx.x],
                data_y_s[threadIdx.y][threadIdx.x],
                data_x_s[threadIdx.y+1][threadIdx.x],
                data_y_s[threadIdx.y+1][threadIdx.x]);

    }
    __syncthreads();

    if (threadIdx.y==0)//the 1st half warp do this.
    {
        float L = 0;
        int i;
        for (i=0; i<=prevP; i++)
        {
            L+=segmentLen[i][threadIdx.x];
        }

        L /= (P + 1);
        int curSegment = 0;
        int nextP = 1;
        float prevSegmentsLen = 0;
        float d;
        for (i = 0; i < P; i++)
        {
            while (segmentLen[curSegment][threadIdx.x] < L * (i + 1) - prevSegmentsLen)
            {
                prevSegmentsLen += segmentLen[curSegment][threadIdx.x];
                curSegment++;
                nextP = curSegment+1;
            }
            d = L * (i + 1) - prevSegmentsLen;

            newEdgePoints_x[i][threadIdx.x]=data_x_s[nextP-1][threadIdx.x] +
                                            (data_x_s[nextP][threadIdx.x] - data_x_s[nextP-1][threadIdx.x]) * d / segmentLen[curSegment][threadIdx.x];

            newEdgePoints_y[i][threadIdx.x]=data_y_s[nextP-1][threadIdx.x] +
                                            (data_y_s[nextP][threadIdx.x] - data_y_s[nextP-1][threadIdx.x]) * d / segmentLen[curSegment][threadIdx.x];
        }
    }

    __syncthreads();

    //write back to global mem, each line of threads handle one column.
    if (tIndex<numEdges&&threadIdx.y<P)
    {
        data_x_d[(threadIdx.y+2)*numEdgesEx+tIndex]=newEdgePoints_x[threadIdx.y][threadIdx.x];
        data_y_d[(threadIdx.y+2)*numEdgesEx+tIndex]=newEdgePoints_y[threadIdx.y][threadIdx.x];
    }
}



//********************************************************************
//the first cycle of bundle.
void firstCycle()
{
    // set parameters for the next cycle.
    int prevP=P;

    dim3 gridDim_amp(1,(480+63)/64);//the dimention of addMidPoints.
    dim3 blockDim_amp(1,64);

    dim3 gridDim_csf(1,480/16);//the dimention of computeSpringForce.
    dim3 blockDim_csf(16,2);

    dim3 gridDim_cef(480/16,480/16);//the dimention of computeElectroForce.
    dim3 blockDim_cef(16,16);

    dim3 gridDim_sef(1,480/16);//the dimention of sumElectroForce.
    dim3 blockDim_sef(16,16);

    dim3 gridDim_us(1,(480+63)/64);//the dimention of updateSubdivs.
    dim3 blockDim_us(1,64);

    addMidPoints<<<gridDim_amp,blockDim_amp>>>(data_x_d,data_y_d);
    cudaThreadSynchronize();


    for(int step=0; step<I; step++)
    {
        computeSpringForce<<<gridDim_csf,blockDim_csf>>>(S,P,edgeLengths_d,data_x_d,data_y_d,Fsi_x_d,Fsi_y_d);
        cudaThreadSynchronize();

        computeElectroForce<<<gridDim_cef,blockDim_cef>>>(P,S,compatibility_d,data_x_d+2*numEdgesEx,data_y_d+2*numEdgesEx,Fei_x_d_0,Fei_y_d_0,edgeValues_d,edgeValueMax,edgeValueMin);
        cudaThreadSynchronize();

        sumElectroForce<<<gridDim_sef,blockDim_sef>>>(P,Fei_x_d_0,Fei_y_d_0,Fei_x_d,Fei_y_d);
        cudaThreadSynchronize();

        updateSubdivs<<<gridDim_us,blockDim_us>>>(P,data_x_d+2*numEdgesEx,data_y_d+2*numEdgesEx,Fsi_x_d,Fsi_y_d,Fei_x_d,Fei_y_d);
        cudaThreadSynchronize();
    }

}



//********************************************************************
//the cycle of bundle other than the first one.
void nextCycle()
{
    // set parameters for the next cycle.
    int prevP=P;

    Pfloat *= Macro_subdivisionPointsCycleIncreaseRate;
    P = (int)(Pfloat+0.5);
    S *= (1.0 - Macro_stepDampingFactor);
    I = (I * 2) / 3;

    int dimy=(P+2)>>1;
    dimy<<=1;

    dim3 gridDim_asp(1,480/16);
    dim3 blockDim_asp(16,dimy);

    dim3 gridDim_csf(1,480/16);
    dim3 blockDim_csf(16,dimy);

    dim3 gridDim_cef(480/16,480/16);
    dim3 blockDim_cef(16,16);

    dim3 gridDim_sef(1,480/16);
    dim3 blockDim_sef(16,16);

    dim3 gridDim_us(1,(480*P+63)/64);
    dim3 blockDim_us(1,64);


    addSubdivisionPoints<<<gridDim_asp,blockDim_asp>>>(P,prevP,data_x_d,data_y_d);
    cudaThreadSynchronize();

    for(int step=0; step<I; step++)
    {
        computeSpringForce<<<gridDim_csf,blockDim_csf>>>(S,P,edgeLengths_d,data_x_d,data_y_d,Fsi_x_d,Fsi_y_d);
        cudaThreadSynchronize();

        computeElectroForce<<<gridDim_cef,blockDim_cef>>>(P,S,compatibility_d,data_x_d+2*numEdgesEx,data_y_d+2*numEdgesEx,Fei_x_d_0,Fei_y_d_0,edgeValues_d,edgeValueMax,edgeValueMin);
        cudaThreadSynchronize();

        sumElectroForce<<<gridDim_sef,blockDim_sef>>>(P,Fei_x_d_0,Fei_y_d_0,Fei_x_d,Fei_y_d);
        cudaThreadSynchronize();

        updateSubdivs<<<gridDim_us,blockDim_us>>>(P,data_x_d+2*numEdgesEx,data_y_d+2*numEdgesEx,Fsi_x_d,Fsi_y_d,Fei_x_d,Fei_y_d);
        cudaThreadSynchronize();
    }
}



//********************************************************************
//bundle.
void bundle()
{

    //allocate mem.
    size_t size=numEdgesEx*(numSubdivs+1)*sizeof(float);//spring force.
    cudaMalloc((void**)&Fsi_x_d,size);
    cudaMalloc((void**)&Fsi_y_d,size);

    size=numEdgesEx*(numSubdivs)*sizeof(float);//electrostatic force.
    cudaMalloc((void**)&Fei_x_d,size);
    cudaMalloc((void**)&Fei_y_d,size);

    float *Fei_x_h=(float *)malloc(size);
    float *Fei_y_h=(float *)malloc(size);

    size=numEdgesEx*numEdgesEx*(numSubdivs)*sizeof(float);//intermediate result of electrostatic force.
    cudaMalloc((void**)&Fei_x_d_0,size);
    cudaMalloc((void**)&Fei_y_d_0,size);

    // iterative refinement scheme.
    firstCycle();
    for (int cycle= 1;cycle < Macro_numCycles; cycle++)
    {
        nextCycle();
    }

    size=numEdgesEx*(numSubdivs)*sizeof(float);
    cudaMemcpy(data_x_h+2*numEdgesEx,data_x_d+2*numEdgesEx,size,cudaMemcpyDeviceToHost);
    cudaMemcpy(data_y_h+2*numEdgesEx,data_y_d+2*numEdgesEx,size,cudaMemcpyDeviceToHost);
}


//********************************************************************
//read the coordinate of the points.
void initPoints()
{
    FILE *f;
    f=fopen("D://cord_cord.txt","r");//FIXME
    float x=0,y=0;
    for(int i=0; i<numPoints; i++)
    {
        fscanf(f,"%f,",&x);
        fscanf(f,"%f",&y);
        endPoints[i].x=x;
        endPoints[i].y=y;
    }
    fclose(f);
}


//********************************************************************
//compute the lengths of the edges.
void initValues()
{
    float evMin = pow(10.0,10.0), evMax = -pow(10.0,10.0);
    FILE *f;
    f=fopen("D://mig2005.txt","r");//FIXME
    int i,j;
    int k;
    int value;
    float length;
    char v[256];

    size_t size;
    size=numEdgesEx*(numSubdivs+2)*sizeof(float);
    data_x_h=(float *)malloc(size);
    data_y_h=(float *)malloc(size);


    size=numEdgesEx*sizeof(float);
    edgeLengths_h=(float*)malloc(size);
    edgeValues_h=(float*)malloc(size);

    fscanf(f,"%s",v);
    fgetc(f);
    for(i=1,k=0; i<numPoints; i++)
    {
        for(j=0; j<i-1; j++,k++)
        {
            data_x_h[k]=endPoints[j].x;
            data_y_h[k]=endPoints[j].y;

            data_x_h[numEdgesEx+k]=endPoints[i].x;
            data_y_h[numEdgesEx+k]=endPoints[i].y;

            length = computeDistance(endPoints[i].x,endPoints[i].y,endPoints[j].x,endPoints[j].y);
            if (abs(length) < EPS) length = 0.0;
            edgeLengths_h[k]=length;

            fscanf(f,"%d",&value);
            fgetc(f);
            edgeValues_h[k]=value;
            if(value>evMax)
            {
                evMax=value;
            }
            if(value<evMin)
            {
                evMin=value;
            }
        }
        data_x_h[k]=endPoints[j].x;
        data_y_h[k]=endPoints[j].y;

        data_x_h[numEdgesEx+k]=endPoints[i].x;
        data_y_h[numEdgesEx+k]=endPoints[i].y;

        length = computeDistance(endPoints[i].x,endPoints[i].y,endPoints[j].x,endPoints[j].y);
        if (abs(length) < EPS) length = 0.0;
        edgeLengths_h[k]=length;

        fscanf(f,"%d",&value);
        fgetc(f);
        fscanf(f,"%s",v);
        fgetc(f);
        edgeValues_h[k]=value;
        if(value>evMax)
        {
            evMax=value;
        }
        if(value<evMin)
        {
            evMin=value;
        }
        k++;
    }
    fclose(f);
    if (Macro_edgeValueAffectsAttraction)
    {
        edgeValueMax = evMax;
        edgeValueMin = evMin;
    }
}


//********************************************************************
//init all the required stuff.
void init()
{
    initPoints();
    initValues();

    I = Macro_I;
    P = Macro_P;
    Pfloat = P;
    S = Macro_S;

    //CPU->GPU
    size_t size=numEdgesEx*(numSubdivs+2)*sizeof(float);
    size_t cp_size=numEdgesEx*2*sizeof(float);
    cudaMalloc((void**)&data_x_d,size);
    cudaMemcpy(data_x_d,data_x_h,cp_size,cudaMemcpyHostToDevice);

    cudaMalloc((void**)&data_y_d,size);
    cudaMemcpy(data_y_d,data_y_h,cp_size,cudaMemcpyHostToDevice);

    //edgeLengths,edgeValues.
    size=numEdgesEx*sizeof(float);
    cudaMalloc((void**)&edgeLengths_d,size);
    cudaMemcpy(edgeLengths_d,edgeLengths_h,size,cudaMemcpyHostToDevice);

    cudaMalloc((void**)&edgeValues_d,size);
    cudaMemcpy(edgeValues_d,edgeValues_h,size,cudaMemcpyHostToDevice);


}


//********************************************************************
//output the result matrix to a file, one can use it to draw a beautiful graph. 
void writeGraph()
{
    FILE *f;
    f=fopen("D://CUDA_graph_bundler.txt","w");//FIXME
    int i;
    float p_x,p_y;
    for(int pe=0; pe<numEdges; pe++)
    {
        p_x = data_x_h[numEdgesEx+pe];
        p_y = data_y_h[numEdgesEx+pe];
        fprintf(f,"%f,%f,",p_x,p_y);
        for(i=0; i<P; i++)
        {
            p_x=data_x_h[(i+2)*numEdgesEx+pe];
            p_y=data_y_h[(i+2)*numEdgesEx+pe];
            fprintf(f,"%f,%f,",p_x,p_y);
        }
        p_x = data_x_h[pe];
        p_y = data_y_h[pe];
        fprintf(f,"%f,%f",p_x,p_y);
        fprintf(f,"\n");
    }
    fclose(f);
}


//********************************************************************
int main()
{
	float elapsedTime;
	cudaEvent_t start,stop;
    
	//input.
	init();
	
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start,0);//record start time.
	
	//bundle.
    calcEdgeCompatibilityMeasures();
    bundle();

	cudaEventRecord(stop,0);//record stop time.
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime,start,stop);
	cudaEventDestroy(start);
	cudaEventDestroy(stop);
	printf("time: %f ms\n",elapsedTime);//output used time.

	//ouput result.
	writeGraph();
	
    return 0;
}



